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Biophysical impacts of Earth greening largely 
controlled by aerodynamic resistance 
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Satellite observations show widespread increasing trends of leaf area index (LAI), known as the Earth greening. 
However, the biophysical impacts of this greening on land surface temperature (LST) remain unclear. Here, we 
quantify the biophysical impacts of Earth greening on LST from 2000 to 2014 and disentangle the contributions 
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of different factors using a physically based attribution model. We find that 93% of the global vegetated area 
shows negative sensitivity of LST to LAI increase at the annual scale, especially for semiarid woody vegetation. 
Further considering the LAI trends (P < 0.1), 30% of the global vegetated area is cooled by these trends and 5% is 
warmed. Aerodynamic resistance is the dominant factor in controlling Earth greening’s biophysical impacts: The 
increase in LAI produces a decrease in aerodynamic resistance, thereby favoring increased turbulent heat transfer 


between the land and the atmosphere, especially latent heat flux. 


INTRODUCTION 

Vegetation is a key regulator of land-atmosphere exchanges of heat, 
mass, and momentum (1, 2). Satellite records of leaf area index 
(LAI), a measure of the aboveground vegetation abundance, indi- 
cate widespread increasing trends since the 1980s (3), due to warm- 
ing (4, 5), CO: fertilization (6), and land management (7, 8). The 
Earth greening directly affects the surface energy budget through 
influencing various biophysical factors: (i) albedo (a), which con- 
trols the fraction of solar radiation absorbed by the surface; (ii) 
aerodynamic resistance (ra), which characterizes the efficiency of 
turbulent transfer of heat and water vapor; (iii) surface resistance 
(rs), which is the additional resistance to water vapor transport 
through the soil and the pores on leaves; and (iv) emissivity (£), 
which is the effectiveness of a surface in emitting and absorbing 
longwave radiation. Changes in these biophysical factors can strongly 
affect the radiometric land surface temperature (LST). 

Despite their importance, the magnitude of these biophysical 
changes associated with Earth greening and their impacts on LST 
remain poorly understood (2, 9-12). First, it is unclear which bio- 
physical factor dominates the effects of Earth greening on LST. A 
global modeling study suggested that changes in aerodynamic resist- 
ance caused by the increasing LAI play a negligible role (9). On the 
other hand, several other studies demonstrate that aerodynamic 
resistance is the most critical biophysical factor that controls the 
LST response to vegetation changes (13-15). Quantifying the bio- 
physical impacts of Earth greening on LST and addressing the dom- 
inant biophysical factor frame the scope of our study. 

Second, the biophysical impacts of Earth greening on LST can- 
not be simply quantified by the observed changes in LST from sat- 
ellites. The observed changes in LST are compound effects of Earth 
greening and large-scale climate change (fig. $1)—here, the Earth 
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greening means an increasing trend in LAI, and the large-scale 
climate change refers to changing atmospheric conditions such as 
rising air temperature. Consequently, the sensitivity of LST to LAI 
inferred from statistical methods, e.g., multiple linear regression 
(10), may suffer causality issues [see, e.g., a recent debate on whether 
greening is the cause or effect of warming in the boreal zone (10-12)]. 
In this study, we do not aim to quantify the interplay between Earth 
greening and rising air temperature. Instead, we aim to quantify the 
impacts of Earth greening on LST independent of those from large- 
scale climate change. 

To do so, we use a physically based attribution method, the two- 
resistance mechanism (TRM) method (14, 16, 17). It quantifies the 
sensitivities of LST to Earth greening through biophysical pathways 
[O, ra, fs €, and ground heat flux (G)] and thus the contributions of 
these biophysical changes to changes in LST. Unlike previous studies 
that compare the effects of changes in sensible and latent heat fluxes 
(9, 10, 18, 19), we take one step further by considering changes in 
aerodynamic resistance (mostly related to surface roughness) and 
surface resistance (mostly related to soil moisture and vegetation 
characteristics), which circumvents the strong correlation between 
sensible and latent heat fluxes through aerodynamic resistance 
(14, 16, 17). 

Changes in biophysical factors are caused not only by the Earth 
greening but also by changing atmospheric conditions. For exam- 
ple, changes in surface resistance could be a result of changes in 
LAI, air temperature, and specific humidity (20). Hence, to quantify 
the biophysical impacts of Earth greening on LST, we need to further 
consider the sensitivities of biophysical factors to changes in LAI 
(e.g. T: The final expression of the greening-induced LST change 
is diagnosed as follows 


bio, LAI _ EY Z 
AT? = AAMs 
(ae (atar) + (Gee) arar) * (See) Gane 
da / \OLAI dra / \OLAI ar, / \OLAI 
ALAL at 
+(3) (athr) + (32) (ate) 
de /\ OLAI aG /\ oLAI (1) 
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where ALAIsa is the satellite-derived long-term trend in LAI (fig. S2A) 


multiplied by the length of study period (15 years); arn] is the sensitivity 
of LST to LAI through biophysical pathways; a os ae os and = 


or x ; a Əra rs oe 
are sensitivities of LST to biophysical factors; and 7, SLAP sap LAP 


and oe are sensitivities of biophysical factors to LAI (Materials 
and Methods). 


RESULTS 
Sensitivity of LST to LAl through biophysical factors 
The sensitivity of LST to LAI diagnosed f from Community Land Model 


; 7 
due to the Earth greening globally (Fig. 1A)—93% of the global 


vegetated area shows T < 0 with an average of —0.36 + 0.22 K m” m” 


(mean + 1 SD, where sD indicates spatial variability). We find that the 
2 


ae 
than those i i high- latitude (—0.34 K m? m`’) and tropical regions 


(-0.29 K m? m`’) (Table 1). We further investigated 2 JAI by classi- 
fying the global vegetation into forests, other woody vegetation 
(OWV; un as shrubs and oo, grasslands, and croplands 
(fig. S2B). 3 ILAT 
(-0.43 K m? m~?), followed by Pesak (-0.36 K m? m~?), but 
weak in forests (—0.23 K m? m”’) (Table 1). Our analysis shows that 
OWV mostly grows in relatively dry regions where a change in LAI 
can substantially alter the efficiency of convective heat transfer be- 
tween the land and atmosphere through changing the aerodynamic 
resistance (fig. S3B). A stronger ty in dry regions: is evidenced 
by the fact that the apami o 


LAT 


T 
the saturation effect—the same change in LAI would lead to a 


smaller change in biophysical factors (a, ra and rs) as well as LST 
where LAI i is high (fig. S4, G to I) (15, 21, 322). Asa result, the mag- 
nitude of 22 in regions with LAI <1 m? m” (-0.45 K m? m°) is 


ƏLAI 


120°E 


180° B 180' 
90°N 


much larger than regions with LAI >4 m? m”? (—0.09 K m? m~’) 
(Table 1). 

Multiple lines of evidence show that our 2E a aa 7 diagnosed from 
CLM5 outputs using the TRM method is robust. First, 347 - (Fig. 1A) 
is almost identical to its counterpart (i.e., Ah = = TAN mae ae Materi- 
als and Methods) approximated using CLM5 sensitivity experiments 
by perturbing LAI climatology (Fig. 1D). This directly demonstrates 
that the TRM method can successfully capture the nonlinear re- 
sponse of LST to LAI simulated by CLM5. Second, given that there 
have been studies reporting CLM’s biases in modeling sensible and 
latent heat fluxes (23, 24), we use the Modern-Era Retrospective 
analysis for Research and Applications Version 2 (MERRA-2) re- 
analysis data and seven fully coupled CMIP5 (the fifth phase of the 
Coupled Model Intercomparison Project) model outputs to com- 
pute the sensitivities of LST to biophysical factors (e.g., 25, and we 
find that the results agree well with those diagnosed from CLM5 
outputs (Fig. 1, B and C, and fig. S5). Third, we also use these LST 
sensitivities to goar te factors computed from MERRA-2 and 

(fig. S6). The coefficient of variation (CV) 


ƏLAI 


of thes 


au! 


T 
Values of CV for o pathway are 44% (albedo), 10% (aerodynam- 


ic resistance), and 22% (surface resistance). The relatively high CV 
in the albedo pathway turns out to be mapen because the ab- 


a 
than those from aerodynamic and surface resistances (fig. S6). We 


stress that different inputs for TRM agree in terms of their signs for 
all biophysical pathways, including the albedo pathway (fig. S6). 


Earth greening cools LST 

The change of LST due, to Earth greening (i.e, ATC’) is the 
product of sensitivity =-7 and the long-term change in LAI (ie., 
ALAL) (Eq. 1). Our results show that AT?! is -0.056 + 0.046 K 
over regions with statistically significant LAI trends, and the area of 
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Fig. 1. Biophysical effects of the Earth greening on LST. (A to C) Sensitivity of T, to LAI 


CLM5 outputs, and (C) CMIP5 multi- model ensemble mean (MMEM) and CLM5 outputs. (D) Sensitivity of T; to LAI 
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LAI-induced cooling (30%) is six times of that of LAI-induced 
warming (5%) (Fig. 1F). The magnitude of AT?! estimated here 
is comparable to that reported by a previous study using coupled 
land-atmosphere simulations (9). On the basis of this agreement, 
but acknowledging that the previous study used a different model, 
we conjecture that the Earth greening affects LST mostly through 
biophysical effects with small contributions from atmospheric feed- 
backs. We translate AT?!°"“! from 2000 to 2014 to changes in energy 
(2.97 x 10°" J), which is more than five times the world’s total primary 
energy supply in 2015 (5.71 x 10% J) (25). From this perspective, the 
Earth greening-induced cooling effect is also much stronger than the 
warming caused by land cover change from 2000 to 2015 (1.21 x 10”’J) 
(19), given that one-third of the global vegetated area shows statis- 
tically significant trends in LAI (7). m 

Because the change in LST is the product of the sensitivity, 5z7p 
and the change in LAI, ALAT,a (Eq. 1), high sensitivity of LST to 
LAI does not Deny imply a large change in LST, and vice versa. 
For example, 2=~ Sàr is large in semiarid regions and/or highlands, such 


io 


as in the western United States, western China, Australia, South 
Africa, and Argentina (Fig. 1A). However, these regions have small 
ALAlsat so that AT ie | is suppressed (Fig. 1F „and fig. S2A). 


. JLAI 
large (Fig. 1A), yet ATOLA is negligible because ALAIsat is small 


(Fig. 1F and fig. S2A). In contrast, in subtropical to temperate re- 
gions of the Northern Hemisphere (eastern China, India, band Eastern 
Europe), AT?! is large (Fig. 1F) due to the strong 2 
and ALATgat (fig. S2A). 


aLAl 


The dominant role of aerodynamic resistance 

We quantify the contributions from different biophysical factors to 
JAT Ta Ts €, and G (Eq. 1). At the annual scale, ra plays 
a dominant role in regulating the biophysical impacts of Earth green- 
ing on LST, while the impacts from e and G are two to three orders of 
magnitude smaller than the other factors and thus negligible (26, 27). 
In CLM5, the r, pathway dominates 82% of the global vegetated area, 
followed by the rs pathway (14.7%) and the a pathway (3.6%) (Fig. 2A). 


Table 1. Mean and SD of the biophysical sensitivity of LST to LAI across bioclimatic regimes. Mean + 1 SD, where SD indicates the spatial variability. OWV, 


other woody vegetation. 


Forests owv Grasslands Croplands All vegetation 
Global —0.23 + 0.12 —0.45 +0.32 —0.36 + 0.23 —0.43 + 0.17 —0.36 + 0.22 
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By annual total precipitation (ATP) 
ATP < 900n mm -0. 29+ 0.09 -0.47 +0. 36 -0 35 £0.30 -0.43 +0.17 -0.40 +0.26 
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Fig. 2. Dominant surface biophysical factors in regulating =- a 


at 


yellow, and green represent the dominance of a, ra, and rs, respectively. The inset shows the areal fraction of dominant aes by biome type. FO, forests; OWV, other 


woody vegetation; GR, grasslands; CR, croplands; All, all vegetation. (B) Attribution o 
diamonds indicate the mean. 
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In terms of the areal fraction, 92% of the forests and 89% of OWV 
show that the rą pathway plays a more prominent role than a and rs, 
while the dominant fraction of rą is reduced to 73% in grasslands 
and 68% in croplands (Fig. 2A). The increase in LAI results in a re- 
duced aerodynamic resistance with a global cooling effect of —0.34 + 
0.22 K m? m”? (Fig. 2B), which, on average, contributes to 77% of 
the variance of the total biophysical impacts on LST. When the land 
is hotter and/or wetter than the atmosphere, a reduced aerodynamic 
resistance can enhance the sensible and/or latent heat transfer from 
the land surface to the atmosphere and hence cools the surface 
(1, 28, 29). It can also lead to a land surface warming with a reversed 
heat and water vapor transfer from the atmosphere to the land sur- 
face when the land is cooler than the atmosphere, which is reflected 
in some northern high-latitude regions (fig. S4B). The dominant 
role of r, is robust across results computed from a hybrid of CLM5, 
MERRA-2, and CMIP5 (fig. S6 and table S2) and is consistent 
with another independent research using a different model [i.e., 
MPI-ESM (Max Planck Institute Earth System Model)] (15). 

Compared to the rą pathway, the r, pathway plays a secondary 
role in cooling LST (—0.06 + 0.09 K m? m7”) (Fig. 2B). In terms of 
the areal fraction, the r, pathway is most significant in 29% of crop- 
lands and 19% of grasslands, but least in forests and OWV (8%) 
(Fig. 2A). This is because although rs is more sensitive to changes in 
LAI than r (fig. S4, H and I), LST is much more sensitive to changes 
in r, (i.e. = > =) in CLM5 (fig. S4, E and F) as well as in MERRA-2 
and CMIP5 (fig. S5, K and L). We stress that a minor role of rs does 
not imply a negligible role of latent heat flux because a change in ra 
can also directly cause a change in latent heat flux (Eq. 4) (14-17, 30). 
Using CLM5’s control and sensitivity runs, we compute the sensi- 
tivities of sensible heat flux and latent heat flux to LAI. We find that 
globally, the sensitivity of latent heat flux (6.9 + 4.9 W m~ per unit 
LAI) is larger than that of sensible heat flux (-3.5 + 4.0 W m”? per 
unit LAI) (fig. S7). This implies that increases in LAI generally cause 
the latent heat flux to increase more strongly than the sensible heat 
flux, which actually shows reductions in most places. The total tur- 
bulent flux or available energy (i.e., the sum of sensible and latent 
heat fluxes) increases with increasing LAI. 

The a pathway dominates only in a few places in the Arctic and 
the Sahel (Fig. 2A). In general, the a pathway leads to a minor 
warming effect of 0.04 + 0.07 K m° m” (Fig. 2B). The finding of 
albedo-related warming, especially in northern high altitudes, is 
consistent with previous studies (2, 19, 31, 32), but its magnitude is 
much smaller than the turbulent cooling effects from ra and rs 
(Fig. 2). Except for evergreen forests, the increase in LAI reinforces 
the fraction of absorbed photosynthetically active radiation (fPAR) 
thus reducing a (33), leading to the warming of LST (fig. S4, A and 
G). However, the increase in LAI in evergreen forests has a negli- 
gible o-related impact due to the fact that fPAR is almost saturated 
in evergreen forests (34, 35) and that the increase in LAI enhances 
the reflectance in near-infrared radiation, which elevates a (fig. S4, 
A and G) (21, 36, 37). 


DISCUSSION 

Causality in the attribution analysis 

One highlight of our study is that it treats changes in the LAI as the 
forcing and changes in the LST as the response, thereby avoiding 
the potential causality issues involved in the coevolution between 
LST and LAI. To illustrate this, we use a statistical method follow- 
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ing (1 0) to calculate the sensitivity of LST to LAI. The new sensitiv- 
ity, w, is calculated using multiple linear regression of LST anom- 
alies against LAI, precipitation, and incoming shortwave radiation 
anomalies from CLM5 (Fig. 1E). This new sensitivity can almost 
reproduce the results of (10). However, this new sensitivity is about 
10 times larger than the sensitivity diagnosed by the physically 
based attribution model (=~), and the two have opposite signs 
in, the northern high latitudes (Fig. 1, A and E). Since both ar and 
zaj are computed from CLM5 outputs, they are expected to agree 
with that evaluated from CLMS sensitivity runs with prescribed LAI 


perturbations (An. Fig. 1D and table S1). However, only the former 


is consistent with this expectation. 

The difference between the two methods is likely caused by the 
air temperature. In our CLM5 simulations and the TRM method, 
air temperature is treated as a forcing for the land surface. On the 
other hand, the regression method does not consider the effects of 
anomalies of air temperature on anomalies of LST independent 
from the anomalies of precipitation and incoming shortwave radia- 
tion (10). As a result, the air temperature effects may be spuriously 
attributed to the albedo effects in the northern high latitudes given 
the strong feedbacks in these regions that link the rising tempera- 
ture with albedo changes. To support this, we analytically computed 
the sensitivity of LST to air temperature (i.e., 219 using the TRM 
method, which shows strong positive values over the northern 
high-latitude regions (fig. S8A). This implies that the LST and air 
temperature in these regions are strongly and positively correlated, 
and the rising air temperature inevitably leads to increased LST. The 
rising air temperature in these regions is more likely caused by 
the increasing atmospheric CO; concentration with a minor con- 
tribution from the Earth greening (3, 38). We conducted another 
multiple linear regression with air temperature anomalies as an 
additional independent variable. The resulting =-| becomes nega- 
tive in most places albeit with a different magnitude compared to 
the TRM method (fig. S8B), which suggests that results from the 
regression method are highly dependent on the selected “indepen- 
dent” forcing variables. 


Independence of the attributable variables 

While the simulated cooling effect from the Earth greening in our 
study is consistent with another modeling study (9), our finding of 
the dominant role of aerodynamic resistance is in contrast with 
their conclusion that it plays a negligible role. The main reasons for 
this difference are their use of 2-m air temperature instead of LST 
and, equally importantly, the independence assumption made by 
any attribution methods using first-order Taylor series expansion. 
By neglecting higher-order and cross-order terms, these attribution 
methods assume that the attributing factors (such as a, ra and rs in 
our study) are independent of each other (16). In (9), the authors 
attributed changes in 2-m air temperature to changes in albedo, 
aerodynamic resistance, latent heat flux, atmospheric shortwave 
transmissivity, near-surface air emissivity, and atmospheric circula- 
tion. Hence, they implicitly assumed that latent heat flux is inde- 
pendent of aerodynamic resistance. However, the dependence of 
latent heat flux on aerodynamic resistance is clear from the well- 
accepted parameterization for latent heat flux (Eq. 4) and has been 
demonstrated empirically using eddy covariance observations from 
AmeriFlux (16). Their assumption of independence between latent 
heat flux and aerodynamic resistance may therefore alter the attri- 
bution results (17). 
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With regard to the stronger sensitivity of LST to changes in ra 
than in rs this is not inconsistent with a previous theoretical study 
(30) showing that evaporation efficiency is the most important 
parameter in controlling LST dynamics. First, their estimate of 
evaporation efficiency is for daytime and relatively moist condi- 
tions (30), while our study includes both daytime and nighttime 
and dry conditions at the annual scale. Second, their evaporation 
efficiency is actually dependent on both aerodynamic resistance 
and surface resistance due to the use of a different parameterization 
of latent heat flux (30). Our analysis shows that the Earth greening 
causes large changes in latent heat flux due to the reduced aero- 
dynamic resistance. 


Concluding remarks 

In this study, we evaluate the biophysical impacts of Earth greening 
on LST using an attribution method based on the surface energy 
balance equation. We find that the widespread Earth greening leads 
to a cooling effect on LST across the globe at the annual scale, which 
is predominantly attributed to the decrease in aerodynamic resistance. 
While these small perturbations in LAI tend to alter turbulent pro- 
cesses more than radiative processes globally, radiative processes 
remain critical in a small proportion of regions in the Arctic and 
some sparsely vegetated areas. Last, the TRM method provides a new 
way to diagnose model outputs and can be further used to evaluate 
whether these biospheric impacts of vegetation would be amplified 
or hindered in future climates. If the Earth greening continues, the 
aerodynamic resistance to turbulent transfer will continue to de- 
crease, resulting in stronger instabilities in the atmospheric boundary 
layer. In the meantime, surface resistance will also decrease, possibly 
leading to more water vapor into the atmosphere thus affecting the 
hydrologic cycle. Whether these effects will be detectable by obser- 
vations and whether the Earth greening can affect other climatic 
processes (e.g., extreme events) remain to be investigated. 


MATERIALS AND METHODS 
Model experiments 
We conduct offline land model simulations using the CLM5, which 
is part of the Community Earth System Model Version 2 (CESM2) 
(39), to study the biophysical impacts of Earth greening on the sur- 
face climate. In other words, our study object is the land surface, not 
the coupled land-atmosphere system. Therefore, our CLM5 ex- 
periments do not explicitly include the feedbacks of LAI changes on 
the large-scale climate. Instead, we assume that the ambient at- 
mosphere already considers the impacts of LAI changes and is es- 
sentially the forcing of the land surface. The atmospheric forcing 
is taken from the third phase of the Global Soil Wetness Project [GSWP3 
(http://hydro.iis.u-tokyo.ac.jp/GSWP3/)]. All simulations are conducted 
at 0.47° x 0.63° resolution with a constant atmospheric CO, con- 
centration of 367 parts per million by volume in 2000 to exclude 
effects from vegetation responses to the rising CO, concentration. 
CLM5 explicitly parameterizes the LAI impacts on surface bio- 
physical factors. According to the CLM5 technical note, albedo is 
influenced by LAI through a two-stream approximation where LAI 
affects extinction and scattering coefficients (like the Beer-Lambert 
law) (39). For aerodynamic resistance, the displacement height and 
roughness lengths are both functions of LAI (22, 39). For the sur- 
face resistance, CLM5 uses the Medlyn’s model in which surface 
resistance is affected by LAI through photosynthesis and plant hy- 
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draulics (39-41). Last, emissivity is related to LAI with a negative 
exponential relationship (39). 

To obtain the sensitivities of LST to biophysical factors and the sen- 
sitivities of biophysical factors to LAI (see Eq. 1), we conduct a suite of 
simulations for the period of 2000 to 2014 with prescribed LAI obtained 
from the Moderate Resolution Imaging Spectroradiometer (MODIS): 

1. Spin-up run (1990-1999): We prescribe the MODIS LAI clima- 
tology, i.e., the LAI values change from month to month but have 
no interannual variabilities or trends. The spin-up run ensures that 
the model achieves equilibrium under the current climate. 

2. Control run (2000-2014): All configurations are the same as 
the spin-up run, but the simulation period is from January 2000 
to December 2014. 

3. Sensitivity runs (2000-2014): Compared to the control run, 
we perturb the MODIS LAI climatology by +2%. There are still no 
interannual variabilities or trends for LAI. 

To compare our results to a previous study using multiple linear 
regression method (10), we conduct a historical LAI run using CLM5: 

4. Historical LAI run (2000-2013): We prescribe monthly MODIS 
LAI to CLM5. MODIS LAI is derived from satellite observations, 
which include interannual variabilities. 


MODIS LAI 

The monthly MODIS LAI (0.5° x 0.5°, 2001-2013) at plant func- 
tional type (PFT) level is downloaded from the University Corpora- 
tion for Atmospheric Research (UCAR) data archive for CESM2 
(https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/Ind/ 
clm2/lai_streams/). For each PFT, we calculate the LAI climatology 
and the perturbed LAI climatology at the monthly scale. 

We use Collection 6 Terra and Aqua MODIS LAI products 
(MOD15A2H and MYD15A2H, available at https://e4ftl01.cr.usgs.gov) 
(33) to obtain the linear trend in annual average LAI (2000-2014). 
We check the quality flag for the LAI products to exclude low-quality 
observations contaminated by clouds, aerosols, shadows, snow, and/or 
ice. We apply the same preprocess procedure as documented in a 
previous study (7) and resample the dataset to 0.5° x 0.5° to con- 
form to the resolution requirement of CLM5’s inputs. We calculate 
the trend in annual average LAI using the Mann-Kendall test (R 
package: https://cran.r-project.org/web/packages/zyp/index.html) 
at P < 0.1. The quality of MODIS LAI was comprehensively evaluated 
against ground-based measurements and through intercomparisons 
with other satellite-retrieved LAI products (42, 43).The trend of 
MODIS LAT is also consistent with AVHRR (Advanced Very-High- 
Resolution Radiometer) LAI and other LAI datasets proved by a 
wide range of previous studies (3, 7). 


Diagnosis of the biophysical impacts on LST 

As a result of postindustrial human activities, the atmospheric CO3 
concentration has increased and the land surface has undergone 
substantial changes. These factors collectively contribute to large- 
scale climate change and the observed widespread vegetation green- 
ing (ie., LAI change) (44). Figure S1 shows that both large-scale 
climate change (the gray paths) and LAI change (the blue paths) can 
result in changes in surface biophysical factors (such as albedo, 
aerodynamic resistance, and surface resistance) (28, 45-49), which 
lastly lead to changes in LST (14, 16), denoted as AT?“ (the gray 
and red paths) and AT?4! (the blue and red paths), respectively. 
In addition, large-scale climate change can cause changes in LST 
without altering the biophysical factors (the pink path, denoted as 
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AT*™"™, fig. S1) (46, 47). For example, when the air temperature 
increases, LST would respond even if all the biophysical factors remain 
the same (20). Therefore, only part of the LST change is due to vege- 
tation biophysical impacts AT?OAL The full expression of A TPi>LAI 
is Eq. 1 in the main text. 


Sensitivities of LST to biophysical factors 
Sensitivities of LST to biophysical factors are analytically computed 
from the outputs of CLM5 prescribed with LAI climatology using 
the TRM method. We further examined these sensitivities using re- 
analysis data (MERRA-2) and outputs from seven fully coupled 
CMIP5 models (ACCESS1-0, ACCESS1-3, CCSM4, HadGEM2-ES, 
IPSL-CM5A-LR, IPSL-CM5A-MR, and NorESM1-M). We choose 
these seven CMIP5 models because only they provide all required 
inputs for the TRM method (16). 

We start with the surface energy balance equation, which is ex- 
pressed as 

Ry = Sin(l-a)+eLin-e0T,* = H+LE+G (2) 

where R, is the net radiation; Sin and Lin are the incoming shortwave 
and longwave radiation, respectively; a and e are the albedo and 
emissivity, respectively; H and LE are the sensible and latent heat 
fluxes, respectively; and eoT,’ is the outgoing longwave radiation 
where o is the Stefan-Boltzmann constant and T, is the LST. Further 
connecting H and LE with LST through the aerodynamic resistance 
(ra) and surface resistance (rs) concepts gives 


H = ERT,- T) 6) 
LE = Eini- (4) 


where p is the air density, cp is the specific heat of air at constant 
pressure, Ly is the latent heat of vaporization, T, is the air potential 
temperature, and q, is the air-specific humidity. The use of aerody- 
namic resistance (r,) and surface resistance (rs) gives rise to the name 
of the attribution method [TRM method (14, 16, 17)]. We note that 
the TRM method calculates the bulk r, and r, using Eqs. 3 and 4. 
Substituting Eqs. 3 and 4 into Eq. 2 yields a nonlinear equation 
for Ts, provided that all other variables are given as inputs. This 
equation is further linearized following previous studies (14, 16, 17), 
so that an analytical expression for T; can be obtained, as follows 


do| Ra- G (q3(Ta)—4a) 
Ts = | A a l F Ta (5) 
where 
Ri, = Sa(l-a)+eLin- E0 Tr (6) 
fo 5 fa 
= #1 +¥nen)| ) 
de 
ò = =e 8 
aT la (8) 
fo = P Cp Ào (9) 
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CpP 

Y= 06221, (10) 
1 

11 

4eo T, (11) 


and e* is the saturation vapor pressure and P is the atmospheric pressure 

at the surface. With the analytical expression for T;, we can take 
aT, ƏT; ƏT, OTs aT; 

the partial derivatives——, Or. Or. de? and == to obtain sensitivities 

of LST to biophysical factors, as follows 


OT, n _ AoSin 
ða  l+f (12) 
OT; E dop Lylgi(Ta)-qa) 1 i 
Ora (ratr)? (+f) 
To 5 fa 2 * P Ly(qal Ta) —Ga) 1 
a Bla) his Oo trd CES 
(13) 
OTs _ ho pLy(qi(Ta) —4a) 1 p 
Ors (ratr)? (+f) T 
x 5 To R* p Ly(q3(Ta) —qa) 1 
fe} 2 n 2 
Vtg + 1) (ra + Ts) (+f) 
OT; _ ho 
0G TF (5) 
* p Ly(qa(Ta) —4a) 
IT, Ao ig a ee 
ge Tage ea +f) (16) 


We also calculated the sensitivity of LST to air temperature 


OT, _ OR, ply oq, 
OT, 1+flOT, (ra+rs)ðTa 
R= G- pLy(qi(Ta) —4a) o| Rè- G- pLv(qi(Ta) —4a) 
(ra + Ys) (ta + Ts) +1 
l+f T (+p? oTa 
(17) 
ƏR 3 dho 3 Z 4a _ 062 of _ Pe 
where ar = 74E o Typ I T Ter I = 8 = aji + 


To use the TRM attri- 


late) lie mt Oe yar? at T, = ale 
bution method, the ini a variables are from CLM5 con- 
trol run, which include incoming shortwave radiation, outgoing 
shortwave radiation, incoming longwave radiation, sensible heat 
flux, latent heat flux, ground heat flux, emissivity, surface pressure, 
and air temperature, specific humidity, and air density at the lowest 
atmospheric model level (about 30 m above the land surface) (16). 
We exclude all negative aerodynamic resistance and/or surface re- 
sistance (physically meaningless) inferred from Eqs. 3 and 4. We 
calculate the sensitivities for each year and take the median of these 
sensitivities from 2000 to 2014 for the consequent analyses. 
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Sensitivities of biophysical factors to LAI 

We use outputs from the CLM5 control run and sensitivity runs to 
estimate the sensitivities of biophysical factors to LAI (Eq. 18). For 
example, the albedo sensitivity to LAI is estimated as 


da = Osen — 
oLAI LA Len — 


Qctl 
LAI 


(18) 


where the subscripts “sen” and “ctl” indicate sensitivity and control 
simulations, respectively. Note that a can be replaced with ra rs €, 
and G. As mentioned in the model experiments, there are two sen- 
sitivity runs (i.e., perturb the MODIS LAI climatology by +2%). For 
each year, we calculate the average sensitivities of biophysical fac- 
tors to LAI estimated from the two combinations of sensitivity runs 
and the control run. We take the median sensitivities from 2000 to 
2014 for the consequent analyses. Because of a lack of analytical 
forms of the sensitivities of biophysical factors to LAI, they can only 
be approximated by offline model simulations to exclude the con- 
founding effects from the changes in atmospheric conditions. 


The robustness of the TRM method and the 

diagnosed sensitivities 

To directly examine thę, validity of the TRM method, first, we com- 
pared the diagnosed 2 IAI S o3 
rectly calculated from CLM5 control and sensitivity runs. Similar to 
the biophysical sensitivities to LAI, the 27 can be calculated as 


ALAI 
AT; _ Tsen e Ts,ctl (19) 
ALAI ~ LAT cen — LAT eu 


Second, we cgp pare the LST sensitivities to biophysical factors 
(i.e. sh , z = nd 5) calculated from CLM5 to those (i) from MERRA-2 
(0.625° x Ô. ny 2000- 2017) and (ii) from seven CMIP5 historical runs 
(resolution varies by model, 2000-2005). This can address whether 
the analytically analyzed sensitivities (i.e., oh oh mr, and$ 25 diagnosed 
from offline CLM5 simulations are consistent Mith those diagnosed 
from reanalysis data and fully coupled CMIP5 simulations. MERRA-2 
is the latest version of reanalysis data produced by the NASA’s Global 
Modelling and Assimilation Office (50), which considers the inter- 
actions between lands, atmosphere, and oceans, as well as assimilates 
observational data (50). CMIP5 historical runs are fully coupled 
simulations of recent past that impose changing conditions consist- 
ent with observations, which includes the effects of anthropogenic 
and volcanic influences and solar activities. 

Third, we calculate th a aT 

aT, IT an 


CMIPS, and CLMS. That is, 3%, 27 : 


da dra 
MERRA-2 and CMIP5, and 2%, 2, 


are smn from 


The 


and JAI aa 


es- 


oi sat 


timated purely from CLM5 outputs. 
Fourth, we estimate the sensitivity of LST to LAI using the iden- 
tical multiple linear regression method described in (10) (denoted 


s 
ƏLAI 


ôT; = a+ b- LAI + c - Precipitation + d-d5SWin (20) 


where ôTs, LAI, Precipitation, and 5SW;, are the annual anomalies 
from the CLM5 historical LAI run. We note that the MODIS LAI pro- 
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duct is the input for, ÇLM5’s historical LAI run in satellite phenology 
mode. Naturally, => is the slope b in Eq. 20. Other coefficients a, c, 
and d are not discussed here. Besides, we conducted an additional mul- 
tiple linear regression that includes the interannual variability of air tem- 


JAI 
io 


Contribution of each biophysical pathway to are 


For each grid cell, the contribution by each biophysical pathway is 
defined as the following 


Contribution; = 
[Ai]? 
[fla) ]* + [flra) ]? + (flrs) ]? + (NG) ]? + (Re)? 


where Contribution; is the contribution by a particular biophysical 
Ts ði ; 
pathway, and f(i) = (zZ >) xi ),i = {0 Fa Fs €, G}. 


qPioLAI 


x 100% eh 


The equivalent energy of A 
We note that the equivalent energy of AT? LAT (denoted as Q) is a 
rough estimate that assumes a constant emissivity of unity all over 
places with significant LAI trends. Thus, Q is calculated as 


Q = EnEdex Aix [o (The Bx ATH)” — o(Ti)'}} (22) 


where A is the area of a pixel that varies by latitude; i denotes each 
vegetated pixel with nonzero OY i a. o is the Stefan-Boltzmann 
constant (5.67x10 * W m” K5; Ti is the climatology LST of the 
ith pixel from 2000 to 2014; t is the number of seconds in a year 
(3.1536x10’ s); N is the total years of the study period, which is 15; 
and n is the sequential number of year (from 1 to 15). 


MODIS land cover type product 

The Collection 5.1 MODIS yearly product provides the land cover 
information at 0.05° x 0.05° known as MCD12C1 (51). The overall 
accuracy is around 75%; for example, misclassifications are reported 
between savannas and woody savannas or between cereal croplands 
and grasslands (51-53). We aggregate the International Geosphere- 
Biosphere Programme classification types provided by MCD12C1 
into four broad biome types: forests, OWV, grasslands, and crop- 
lands. Forests consist of evergreen needleleaf forest, evergreen 
broadleaf forest, deciduous needleleaf forest, deciduous broadleaf 
forest, and mixed forest. OWV refers to closed shrublands, open 
shrublands, and woody savannas. Grasslands include savannas and 
grasslands. Croplands consist of croplands and croplands/natural 
vegetation mosaic. There are 12 such global maps, 1 for each year 
from 2001 to 2012. We refine the 12 maps to 1 map by taking the 
mode class of each grid cell. Last, we convert the spatial resolution 
to 0.5° x 0.5° (fig. S2B). 


SUPPLEMENTARY MATERIALS 


Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ 
content/full/6/47/eabb1981/DC1 


REFERENCES AND NOTES 
1. P.J. Sellers, R. E. Dickinson, D. A. Randall, A. K. Betts, F. G. Hall, J. A. Berry, G. J. Collatz, 
A. S. Denning, H. A. Mooney, C. A. Nobre, N. Sato, C. B. Field, A. Henderson-Sellers, 
Modeling the exchanges of energy, water, and carbon between continents 
and the atmosphere. Science 275, 502-509 (1997). 


7 of 9 


ZZOT ‘t0 19q0190 UO SIO‘aNUNTOS' MMM //:sdyy woy popeoumoq 


SCIENCE ADVANCES | RESEARCH ARTICLE 


10. 


11. 


12. 


13. 


14. 


15. 


16. 
17. 
18. 
19. 


20. 


21. 


22. 


23. 


24. 


25. 


Chen etal., Sci. Adv. 2020; 6 : eabb1981 


G. B. Bonan, Forests and climate change: Forcings, feedbacks, and the climate benefits 

of forests. Science 320, 1444-1449 (2008). 

S. Piao, X. Wang, T. Park, C. Chen, X. Lian, Y. He, J. W. Bjerke, A. Chen, P. Ciais, 

H. Tommervik, R. R. Nemani, R. B. Myneni, Characteristics, drivers and feedbacks of global 
greening. Nat. Rev. Earth Environ. 1, 14-27 (2020). 

R. B. Myneni, C. D. Keeling, C. J. Tucker, G. Asrar, R. R. Nemani, Increased plant growth 

in the northern high latitudes from 1981 to 1991. Nature 386, 698-702 (1997). 

L. Xu, R. B. Myneni, F. S. Chapin III, T. V. Callaghan, J. E. Pinzon, C. J. Tucker, Z. Zhu, J. Bi, 

P. Ciais, H. Tommervik, E. S. Euskirchen, B. C. Forbes, S. L. Piao, B. T. Anderson, S. Ganguly, 
R. R. Nemani, S. J. Goetz, P. S. A. Beck, A. G. Bunn, C. Cao, J. C. Stroeve, Temperature 

and vegetation seasonality diminishment over northern lands. Nat. Clim. Chang. 3, 
581-586 (2013). 

Z. Zhu, S. Piao, R. B. Myneni, M. Huang, Z. Zeng, J. G. Canadell, P. Ciais, S. Sitch, 

P. Friedlingstein, A. Arneth, C. Cao, L. Cheng, E. Kato, C. Koven, Y. Li, X. Lian, Y. Liu, R. Liu, 
J. Mao, Y. Pan, S. Peng, J. Peñuelas, B. Poulter, T. A. M. Pugh, B. D. Stocker, N. Viovy, 

X. Wang, Y. Wang, Z. Xiao, H. Yang, S. Zaehle, N. Zeng, Greening of the Earth and its 
drivers. Nat. Clim. Chang. 6, 791-795 (2016). 

C. Chen, T. Park, X. Wang, S. Piao, B. Xu, R. K. Chaturvedi, R. Fuchs, V. Brovkin, P. Ciais, 

R. Fensholt, H. Tommervik, G. Bala, Z. Zhu, R. R. Nemani, R. B. Myneni, China and India 
lead in greening of the world through land-use management. Nat. Sustain. 2, 122-129 
(2019). 

N. Zeng, F. Zhao, G. J. Collatz, E. Kalnay, R. J. Salawitch, T. O. West, L. Guanter, Agricultural 
green revolution as a driver of increasing atmospheric CO; seasonal amplitude. Nature 
515, 394-397 (2014). 

Z. Zeng, S. Piao, L. Z. X. Li, L. Zhou, P. Ciais, T. Wang, Y. Li, X. Lian, E. F. Wood, 

P. Friedlingstein, J. Mao, L. D. Estes, R. B. Myneni, S. Peng, X. Shi, S. I. Seneviratne, Y. Wang, 


Climate mitigation from vegetation biophysical feedbacks during the past three decades. 


Nat. Clim. Chang. 7, 432-436 (2017). 

G. Forzieri, R. Alkama, D. G. Miralles, A. Cescatti, Satellites reveal contrasting responses 
of regional climate to the widespread greening of Earth. Science 356, 1180-1184 (2017). 
Y. Li, Z. Zeng, L. Huang, X. Lian, S. Piao, Comment on “Satellites reveal contrasting 
responses of regional climate to the widespread greening of Earth”. Science 360, 
eaap7950 (2018). 

G. Forzieri, R. Alkama, D. G. Miralles, A. Cescatti, Response to Comment on “Satellites 
reveal contrasting responses of regional climate to the widespread greening of Earth”. 
Science 360, eaap9664 (2018). 

X. Lee, M. L. Goulden, D. Y. Hollinger, A. Barr, T. A. Black, G. Bohrer, R. Bracho, B. Drake, 
A. Goldstein, L. Gu, G. Katul, T. Kolb, B. E. Law, H. Margolis, T. Meyers, R. Monson, 

W. Munger, R. Oren, K. T. Paw U, A. D. Richardson, H. P. Schmid, R. Staebler, S. Wofsy, 

L. Zhao, Observed increase in local cooling effect of deforestation at higher latitudes. 
Nature 479, 384-387 (2011). 

W. Liao, A. J. Rigden, D. Li, Attribution of local temperature response to deforestation. 
J. Geophys. Res. Biogeosci. 123, 1572-1587 (2018). 

J. Winckler, C. H. Reick, R. M. Bright, J. Pongratz, Importance of surface roughness 

for the local biogeophysical effects of deforestation. J. Geophys. Res. Atmos. 124, 
8605-8618 (2019). 

A. J. Rigden, D. Li, Attribution of surface temperature anomalies induced by land use 
and land cover changes. Geophys. Res. Lett. 44, 6814-6822 (2017). 

D. Li, W. Liao, A. J. Rigden, X. Liu, D. Wang, S. Malyshev, E. Shevliakova, Urban heat island: 
Aerodynamics or imperviousness? Sci. Adv. 5, eaau4299 (2019). 

Y. Li, M. Zhao, S. Motesharrei, Q. Mu, E. Kalnay, S. Li, Local cooling and warming effects 
of forests based on satellite observations. Nat. Commun. 6, 6603 (2015). 

G. Duveiller, J. Hooker, A. Cescatti, The mark of vegetation change on Earth's surface 
energy balance. Nat. Commun. 9, 679 (2018). 

P. Wang, D. Li, W. Liao, A. Rigden, W. Wang, Contrasting evaporative responses 

of ecosystems to heatwaves traced to the opposing roles of vapor pressure deficit 

and surface resistance. Water Resour. Res. 55, 4550-4563 (2019). 

D. Huang, Y. Knyazikhin, W. Wang, D. W. Deering, P. Stenberg, N. Shabanov, B. Tan, 

R. B. Myneni, Stochastic transport theory for investigating the three-dimensional canopy 
structure from space measurements. Remote Sens. Environ. 112, 35-50 (2008). 

X. Zeng, A. Wang, Consistent parameterization of roughness length and displacement 
height for sparse and dense canopies in land models. J. Hydrometeorol. 8, 730-737 
2007). 

P. A. Dirmeyer, L. Chen, J. Wu, C.-S. Shin, B. Huang, B. A. Cash, M. G. Bosilovich, 

S. Mahanama, R. D. Koster, J. A. Santanello, M. B. Ek, G. Balsamo, E. Dutra, 

D. M. Lawrence, Verification of land-atmosphere coupling in forecast models, 
reanalyses, and land surface models using flux site observations. J. Hydrometeorol. 19, 
375-392 (2018). 

. N. Williams, Y. Lu, L. M. Kueppers, W. J. Riley, S. C. Biraud, J. E. Bagley, M. S. Torn, 
Land-atmosphere coupling and climate prediction over the U.S. Southern Great Plains. 
J. Geophys. Res. Atmos. 121, 12125-12144 (2016). 

nternational Energy Agency, Key World Energy Statistics 2017 (OECD, 2017). 


20 November 2020 


26. 


27. 


28. 


29. 


30. 


31; 


32. 


33. 


34. 


35. 


36. 


37. 


38. 


39. 


40. 


41. 


42. 


43. 


44. 


45. 
46. 


47. 
48. 


L. Gu, T. Meyers, S. G. Pallardy, P. J. Hanson, B. Yang, M. Heuer, K. P. Hosman, Q. Liu, 

J. S. Riggs, D. Sluss, S. D. Wullschleger, Influences of biomass heat and biochemical 
energy storages on the land surface fluxes and radiative temperature. J. Geophys. Res. 
112, D02107 (2007). 

J.-Y. Juang, G. Katul, M. Siqueira, P. Stoy, K. Novick, Separating the effects of albedo 
from eco-physiological changes on surface temperature along a successional 
chronosequence in the southeastern United States. Geophys. Res. Lett. 34, L21408 (2007). 
J. Kala, M. Decker, J.-F. Exbrayat, A. J. Pitman, C. Carouge, J. P. Evans, G. Abramowitz, 

D. Mocko, Influence of leaf area index prescriptions on simulations of heat, moisture, 
and carbon fluxes. J. Hydrometeorol. 15, 489-503 (2014). 

S. Irmak, D. Mutiibwa, On the dynamics of canopy resistance: Generalized linear 
estimation and relationships with primary micrometeorological variables. Water Resour. 
Res. 46, W08526 (2010). 


S. M. Bateni, D. Entekhabi, Relative efficiency of land surface energy balance components. 


Water Resour. Res. 48, W04510 (2012). 

R. Alkama, A. Cescatti, Biophysical climate impacts of recent changes in global forest 
cover. Science 351, 600-604 (2016). 

R. M. Bright, E. Davin, T. O'Halloran, J. Pongratz, K. Zhao, A. Cescatti, Local temperature 
response to land cover and management change driven by non-radiative processes. 
Nat. Clim. Chang. 7, 296-302 (2017). 

R. B. Myneni, S. Hoffman, Y. Knyazikhin, J. L. Privette, J. Glassy, Y. Tian, Y. Wang, X. Song, 
Y. Zhang, G. R. Smith, A. Lotsch, M. Friedl, J. T. Morisette, P. Votava, R. R. Nemani, 

S. W. Running, Global products of vegetation leaf area and fraction absorbed PAR 

from year one of MODIS data. Remote Sens. Environ. 83, 214-231 (2002). 

R. B. Myneni, D. L. Williams, On the relationship between FAPAR and NDVI. Remote Sens. 
Environ. 49, 200-211 (1994). 

Y. Tian, R. E. Dickinson, L. Zhou, X. Zeng, Y. Dai, R. B. Myneni, Y. Knyazikhin, X. Zhang, 

M. Friedl, H. Yu, W. Wu, M. Shaikh, Comparison of seasonal and spatial variations of leaf 
area index and fraction of absorbed photosynthetically active radiation from Moderate 
Resolution Imaging Spectroradiometer (MODIS) and Common Land Model. J. Geophys. 
Res. 109, D01103 (2004). 

Y. Tian, R. E. Dickinson, L. Zhou, R. B. Myneni, M. Friedl, C. B. Schaaf, M. Carroll, F. Gao, 
Land boundary conditions from MODIS data and consequences for the albedo 

of a climate model. Geophys. Res. Lett. 31, L05504 (2004). 

R. A. McPherson, A review of vegetation—Atmosphere interactions and their influences 
on mesoscale phenomena. Prog. Phys. Geogr. Earth Environ. 31, 261-285 (2007). 

N. L. Bindoff, P. A. Stott, K. M. AchutaRao, M. R. Allen, N. Gillett, D. Gutzler, K. Hansingo, 
G. Hegerl, Y. Hu, S. Jain, I. I. Mokhov, J. Overland, J. Perlwitz, R. Sebbari, X. Zhang, 
Detection and Attribution of Climate Change: from Global to Regional, in Climate Change 
2013: The Physical Science Basis. Contribution of Working Group | to the Fifth Assessment 
Report of the Intergovernmental Panel on Climate Change (IPCC, 2013), T. F. Stocker, D. Qin, 
G-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex, P. M. Midgley, 
Eds. (Cambridge Univ. Press, Cambridge and New York, 2013). 

D. Lawrence, R. Fisher, C. Koven, K. Oleson, S. Swenson, M. Vertenstein, B. Andre, 

G. Bonan, B. Ghimire, L. van Kampenhout, D. Kennedy, E. Kluzek, R. Knox, P. Lawrence, 
F. Li, H. Li, D. Lombardozzi, Y. Lu, J. Perket, W. Riley, W. Sacks, M. Shi, W. Wieder, C. Xu, 

A. Ali, A. Badger, G. Bisht, P. Broxton, M. Brunke, J. Buzan, M. Clark, T. Craig, K. Dahlin, 

B. Drewniak, L. Emmons, J. Fisher, M. Flanner, P. Gentine, J. Lenaerts, S. Levis, Technical 
Description of version 5.0 of the Community Land Model (CLM) (National Center for 
Atmospheric Research (NCAR), 2018). 

B. E. Medlyn, R. A. Duursma, D. Eamus, D. S. Ellsworth, I. C. Prentice, C. V. M. Barton, 

K. Y. Crous, P. de Angelis, M. Freeman, L. Wingate, Reconciling the optimal 

and empirical approaches to modelling stomatal conductance. Glob. Chang. Biol. 17, 
2134-2144 (2011). 

G. B. Bonan, M. Williams, R. A. Fisher, K. W. Oleson, Modeling stomatal conductance 

in the earth system: Linking leaf water-use efficiency and water transport along 

the soil-plant-atmosphere continuum. Geosci. Model Dev. 7, 2193-2222 (2014). 

K. Yan, T. Park, G. Yan, Z. Liu, B. Yang, C. Chen, R. R. Nemani, Y. Knyazikhin, R. B. Myneni, 
Evaluation of MODIS LAI/FPAR Product Collection 6. Part 2: Validation 

and Intercomparison. Remote Sens. 8, 460 (2016). 

K. Yan, T. Park, G. Yan, C. Chen, B. Yang, Z. Liu, R. R. Nemani, Y. Knyazikhin, R. B. Myneni, 
Evaluation of MODIS LAI/FPAR Product Collection 6. Part 1: Consistency and improvements. 
Remote Sens. 8, 359 (2016). 

M. Heimann, M. Reichstein, Terrestrial ecosystem carbon dynamics and climate 
feedbacks. Nature 451, 289-292 (2008). 

W. Brutsaert, Hydrology (Cambridge Univ. Press, 2005). 

R. B. Stull, An Introduction to Boundary Layer Meteorology (Springer Science & Business 
Media, 1988). 

J. R. Garratt, The Atmospheric Boundary Layer (Cambridge Univ. Press, 1994). 

P.G. Jarvis, The interpretation of the variations in leaf water potential and stomatal 
conductance found in canopies in the field. Philos. Trans. R. Soc. Lond. B 273, 593-610 
(1976). 


8 of 9 


ZZOT “PO 18q0190 UO SIO‘aNUNTOS' mMM //:sdyy woy popeoumoq 


SCIENCE ADVANCES | RESEARCH ARTICLE 


49. J.B. Stewart, Modelling surface conductance of pine forest. Agric. For. Meteorol. 43, 19-35 
(1988). 

50. R. Gelaro, W. M. Carty, M. J. Suárez, R. Todling, A. Molod, L. Takacs, C. Randles, A. Darmenov, 
M. G. Bosilovich, R. Reichle, K. Wargan, L. Coy, R. Cullather, C. Draper, S. Akella, V. Buchard, 
A. Conaty, A. da Silva, W. Gu, G.-K. Kim, R. Koster, R. Lucchesi, D. Merkova, J. E. Nielsen, 
G. Partyka, S. Pawson, W. Putman, M. Rienecker, S. D. Schubert, M. Sienkiewicz, B. Zhao, 
The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). 
J. Climate 30, 5419-5454 (2017). 

51. M.A. Friedl, D. K. Mclver, J. C. F. Hodges, X. Y. Zhang, D. Muchoney, A. H. Strahler, 
C. E. Woodcock, S. Gopal, A. Schneider, A. Cooper, A. Baccini, F. Gao, C. Schaaf, Global 
land cover mapping from MODIS: Algorithms and early results. Remote Sens. Environ. 83, 
287-302 (2002). 

52. C. Giri, Z. Zhu, B. Reed, A comparative analysis of the Global Land Cover 2000 and MODIS 
land cover data sets. Remote Sens. Environ. 94, 123-132 (2005). 

53. S. Fritz, L. See, |. McCallum, C. Schill, M. Obersteiner, M. van der Velde, H. Boettcher, 
P. Havlík, F. Achard, Highlighting continued uncertainty in global land cover maps 
for the user community. Environ. Res. Lett. 6, 044005 (2011). 


Acknowledgments 

Funding: D.L. and M.H. were supported by the U.S. Department of Energy, Office of Science, 
as part of research in Multi-Sector Dynamics, Earth and Environmental System Modelling 
Program. P.G. acknowledges funding from the NSF award number 1734156. R.B.M. 


Chen etal., Sci. Adv. 2020; 6 : eabb1981 20 November 2020 


acknowledges funding by the Alexander von Humboldt Foundation. We thank Boston 
University and the National Center for Atmospheric Research (supported primarily by the 
NSF) for providing supercomputing resources. Author contributions: D.L. and C.C. 
conceived the methods and experiments. C.C. carried out model simulations and data 
analyses. C.C. and D.L. wrote the article. All authors contributed comments and critiques on 
drafts. Competing interests: The authors declare that they have no competing interests. 
Data and materials availability: All raw data analyzed in this study are publicly available as 
referenced within the article. CESM2/CLM5 release code is available at www.cesm.ucar.edu/ 
models/cesm2/release_download.html. The trend test package is available at https:// 
cran.r-project.org/web/packages/zyp/index.html. The scripts of the TRM algorithm are 
available at https://github.com/yaoganchenchi/TRM_Science_Advances. Additional data 
related to this paper may be requested from the authors. 


Submitted 6 February 2020 
Accepted 7 October 2020 
Published 20 November 2020 
10.1126/sciadv.abb1981 


Citation: C. Chen, D. Li, Y. Li, S. Piao, X. Wang, M. Huang, P. Gentine, R. R. Nemani, R. B. Myneni, 
Biophysical impacts of Earth greening largely controlled by aerodynamic resistance. Sci. Adv. 6, 
eabb1981 (2020). 


9of9 


TTOT ‘t0 19q0190 UO JI0'avuarnos mMM //:sdyy woy popeoumoq 


Science Advances 


Biophysical impacts of Earth greening largely controlled by aerodynamic 
resistance 


Chi ChenDan LiYue LiShilong PiaoXuhui WangMaoyi HuangPierre GentineRamakrishna R. NemaniRanga B. Myneni 


Sci. Adv., 6 (47), eabb1981. + DOI: 10.1126/sciadv.abb1981 


View the article online 


https:/Awww.science.org/doi/10.1126/sciadv.abb1981 
Permissions 


https://www.science.org/help/reprints-and-permissions 


Use of this article is subject to the Terms of service 


Science Advances (ISSN 2375-2548) is published by the American Association for the Advancement of Science. 1200 New York Avenue 


NW, Washington, DC 20005. The title Science Advances is a registered trademark of AAAS. 
Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim 
to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC). 


ZZOT “PO 19q0190 UO SIO‘a0UNTIOS' MMM //:sdyy woy popeoumoq 


